Molecular-dynamics thermal annealing model 
of laser ablation of silicon 
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A molecular-dynamics thermal annealing model is proposed to investigate the mechanisms in- 
volved in picosecond pulsed laser ablation of crystalline silicon. In accordance with the thermal an- 
nealing model, a detailed description of the microscopic processes which result from the interaction 
of a 308 nm, 10 ps, Gaussian pulse with a Si(100) substrate has been embedded into a molecular- 
dynamics scheme. This was accomplished by explicitly accounting for carrier-phonon scattering and 
carrier diffusion. Below the predicted threshold fluence for ablation, F t n = 0.25 J/cm 2 , a surface 
breathing mode indicates that the solid restores internal equilibrium by the generation of pressure 
waves. Above F t h, our simulations reveal that matter removal is triggered by subsurface superheat- 
ing effects: intense heating of the material leads to the thermal confinement of the laser-deposited 
energy. As a result, the material is overheated up to a temperature corresponding to the critical 
point of silicon and a strong pressure gradient builds up within the absorbing volume. At the same 
time, diffusion of the carriers into the bulk leads to the development of a steep temperature gradient 
beneath the surface. Matter removal is subsequently driven by the relaxation of the pressure gradi- 
ent: large pieces — several atomic layers thick — of molten material are expelled from the surface 
with initial axial velocities of ~ 1000 m/s, their ejection following the nucleation of voids beneath 
the surface. 

PACS numbers: 79.20.Ds, 61.80.Az, 79.20.Ap, 81.15.Fg 
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I. INTRODUCTION 

The use of photons foe, the controlled removal of mat- 
ter, i.e., laser ablation,BEI has given rise to a large num- 
ber of applications for materials processing and growth; 
in microelectronics, examples include surface microma- 
chining. surface cleaning, and the pulsed laser deposition 
(PLD)u of thin films n which allows a wide variety of ma- 
terials to be grown. Bel In spite of this, a complete picture 
of the mechanisms underlying ablation for a broad range 
of pulse durations and wavelengths, and for such diverse 
materials as metals, polymers and semiconductors, is still 
lacking. 

The great complexity of phenomena involved in laser 
ablation, determined by the coupling of the laser pulse 
with the optical, elastic and thermal properties of the 
material, has been an obstacle in understanding. On 
a mesoscopic scale, the processes brought about by ul- 
trashort -^-femtosecond — pulses have seldom been 
addressed. EfflQ On the other haad, evidence of a conmfe 
tition between photomechanicatm and photothermalM2l 
mechanisms has begun to emerge in the light-induced re- 
moval of matter with picosecond and nanosecond pulses. 
In the stress confinement regimejjt] ejection of matter 
proceeds from photomechanical effects: energy is de- 
posited in the material in a time shorter than that needed 
for the generation of acoustic waves, and ablation is in- 
duced by the relaxation of the pressure gradient that 



builds up within the absorbing volume. In the ther- 
mal confinement regime,Q energy is supplied to the ir- 
radiated volume at a rate faster than that at which it 
can be pumped out by thermal diffusion. Consequently, 
the material is rapidly brought near its critical point 
as matter removal is driven by photothermal effects :E£l 
the metastable superheated liquid undergoes a prompt 
transition to a gaseous phase consisting, to a large ex- 
tent, of molten micron-sized droplets. This mechanism 
is known as explosive boiling or phase explosion other 
photothermal processes include vaporization and normal 
boiling (see Ref. |ll| and refer.encf?.s therein). Finally, there 
have been repeated reportsOElilEaE^I of a fourth pro- 
cess known as subsurface superheating (SSSH). In this 
scenario, heating of the material causes evaporation as 
a result of which a steep temperature gradient develops 
beneath the surface. Matter removal then takes place in 
a manner similar to a phase explosion following the re- 
laxation of the pressure build-up in the bulk. It is not 
clear whether this mechanism | is ,| real|as conflicting views 
can be found in the literature. ll3lliHlall3 

While it is generally assumed that the laser sputtering 
of metal surfaces involves phase explosion mechanisms 
(see, e.g., Ref. [n]), both photomechanical and photother- 
mal effects have been shown to be operative in organic 
solids.u In semiconductor materials such as silicon, how- 
ever, understanding of matter removal processes induced 
by picosecond pulses has yet to come. In the present 
study, we give evidence of a SSSH effect in the laser ab- 
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lation of silicon. This is achieved by carrying out detailed 
molecular-dynamics (MD) simulations of the interaction 
of a picosecond laser pulse with a Si(100) substrate, ex- 
plicitly accounting for carrier-phonon scattering and car- 
rier diffusion. _Il is assumed that the thermal annealing 
model (TAM)Ea is valid so that thermal processes only 
are operative. In our computer code, 'MADTAM', we 
use the Stillinger- Weber (SW)lij potential to model the 
atomic interactions; preliminary results from an earlier 
version of our model can be found in Refs. ^ and 21 



Anticipating our results, we find that, below the (pre- 
dicted) threshold fluence for ablation, Fth, the solid re- 
stores internal equilibrium by the generation of pressure 
waves. Above F t h, our simulations reveal that matter re- 
moval is triggered by SSSH effects: intense heating of the 
material leads to the thermal confinement of the laser- 
deposited energy and a steep temperature gradient be- 
neath the surface. Interestingly, the latter is not caused 
by the evaporation of atoms from the outer surface but, 
rather, originates from the process of carrier diffusion into 
the bulk. As a result, a strong pressure gradient builds 
up ~ 30 nm below the surface as the material is over- 
heated up to a temperature corresponding to the critical 
point of silicon. Ablation is subsequently initiated by the 
relaxation of the pressure gradient: large pieces — sev- 
eral atomic layer thick — of molten material are expelled 
from the surface, a consequence of the nucleation of voids 
in the bulk. Before discussing these results, however, we 
give a detailed description of our model. 



II. THE MODEL 

The processing of silicon surfaces with femtosecond 
laser pulses involves ultrafast, non-thermal, phenomena 
which cannot be accounted for by the SW potential. 
For this reason, we remain within the picosecond regime 
where thermal mechanisms are involved; this is discussed 
in Sec. II A. A description of the computational method, 
i.e., molecular dynamics, follows in Sec. II B. Finally, 
the generation of electron- hole (e-h) pairs by photons, as 
well as their relaxation through carrier-phonon scatter- 
ing and carrier diffusion, depicted in Fig. [j], are detailed 
in Sec. ET3toCTE. 



A. Preliminary considerations 

Various computational techniques allow 
the intemrtiniji_ 1| -of_ 1| light .--with matter to be 
simulatedBByyid-BBBm-El Among them, 
MD appears to be a most suitable approach for ad- 
dressing the numerous, complex processes involved 
in the laser ablation and desorption of matter: no 
hypotheses are needed to account for the various phe- 
nomena and atomic motion can be followed, in real time, 
provided that the, interactions between the atoms are 
correctly modeledEilH However, because laser ablation 



involves large length and time scales, and thus a great 
computation al .. effort , ■ .^ temp ts involving classical MD 
are scarce. B^ESEZISEj'ESE^ In order to extend the 
range of MD simulations, Zhigilei et al. proposed a 
breathing-sphere model appropriate to the study of 
matter removal in organic solids.c£l In this work, the 
phase explosion due to overheating and the laser-induced 
pressure build-up were identified as the key processes 
responsible for ablation in the regimes of thermal and 
stress confinements, respectively.0E3 

Another MD Stwdy of laser ablation was performed by 
Herrmann et al.,c3 who examined the interaction of fem- 
tosecond laser pulses with silicon; a truncated version of 
the SW potential was used to account for the rupture of 
interatomic bonds following the absorption of photons. 
The predicted values of the threshold fluence for abla- 
tion thus obtained were at least an order .of. magnitude 
higher than those observed experimentally.EflE3 This dis- 
crepancy is likely due to the inability of the SW potential 
to describe the non-thermal mechanisms involved in the 
interaction of femtosecond pulses with silicon. Indeed, 
it is generally agreecOL3Ej that there exists, for silicon, 
a critical carrier density n c ~ 10 22 cm~ 3 separating two 
distinct regimes: 

(i) With picosecond and nanosecond pulses, carrier 
densities remain below the critical value and the relax- 
ation is thermal, i.e., the carriers relax by transferring 
their kinetic energy to the lattice by the spontaneous 
emission of optical phonons in a characteristic time of 
~ 1 ps. In this context, structural modifications such as 
melting and ablation are observed on a ~ 100 — 1000 ps 
time scale. This model, known as TAM,.is suitable for 
laser pulses down to approximately 10 psJlj 

(ii) With femtosecond pulses, on the other hand, 
the number of free carriers may exceed n c and ultra- 
fast melting-.pijoceeds within ~ 1 ps by the collapse of 
the lattice.E3c3 Because melting occurs before the car- 
riers have time to relax through the emission of optical 
phonons, the relaxation is non-thermal. This model, first 
proposed by Van Vechten-e£ al.,c3 is known as the plasma 
annealing model (PAM)Jl3 

Furthermore, if the light intensity at the surface ex- 
ceeds ~ 10 4 GW/cm 2 , the generation of a high num- 
ber of energetic carriers through a two-photon absorp- 
tion mechanism leads to the complete ionization of the 
material by avalanche ionization and/or field ionization: 
a very dense (n > 10 23 cm _3 ) r and.hot (~ 10 6 K) e-h gas 
is generated near the surf ace. CJ'cfJ One then speaks of a 
solid-to-plasma transition as the material is carried away 
by the expanding plasma. 

Here, because we make use of the SW potential, it is 
assumed that thermal processes only are operative; this 
restricts us to 10 ps or longer pulses. Care is taken to 
avoid carrier densities in excess of n c ; in this context, the 
intensity remains well below ~ 10 4 GW/cm 2 . 



B. Molecular dynamics 
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The MD technique has been discussed extensively in 
the literature (see, for instance, Ref. [52|). Briefly, the set 
of 3iV Newton equations of motion, where N is the num- 
ber of atoms in the system, are solved from the knowl- 
edge of the interatomic forces. If r,;, Vi, and are the 
position, velocity and acceleration of the ith atom, re- 
spectively, the time evolution of the system is obtained 
by integrating, at discrete time steps, the following dif- 
ferential equations: 



dvi 
~dt 

dt 



1 



dUjri, ...,r N ) 
dr 



(1) 



given a set of initial positions and velocities {r°,v?}. In 
the above equations, m, is the mass of the ith atom and 
U the interatomic potential. Thermodynamic quantities 
such as temperature and pressure can then be computed 
in a straightforward manner .til 

All calculations were performed using the program 
groF, a general-purpose MD code for bulk and surfaces 
developed by one of the authors (LJL). The original form 
of the SW potential was-used to compute the forces be- 
tween the silicon atoms£3 The equations of piption were 
solved using the velocity Verlet algorithmic The MD 
time step was set to At — 0.5 fs for all simulations. 



C. Target and laser pulse 

Simulations are carried out for a supercell with approx- 
imate dimensions 3 x 3 x 60 nm 3 (containing a total of 
31680 atoms) representing a small volume of a Si(100) 
target located at the center of the laser pulse. 

In order to recreate the thermal and structural con- 
straints of a macroscopic crystal, the supercell is first re- 
peated in the lateral (x and y) directions using periodic 
boundary conditions (PBC). The thermal constraint is 
imposed by coupling a few monolayers of atoms at the 
bottom of the supercell to a heat reservoir; this is done by 
renormalizing the velocities to an appropriate Maxwell- 
Boltzmann (MB) distribution. These atoms also serve to 
minimize the reflection of pressure waves at the bottom of 
the supercell. This latter phenomenon affects the results 
by, in some cases, adding an artificial contribution to ab- 
lation. Though other techniques were proposed to atten- 
uate such pressure wavesjla the method described above 
gives satisfactory results. Additionally, a few monolay- 
ers of atoms attached to their equilibrium positions are 
placed beneath the heat reservoir in order to mimic a 
semi-infinite crystal. The crystal is assumed to be ini- 
tially perfect. 

The simulation of a 10 ps laser pulse, Gaussian in time 
as well as in the x and y directions and of macroscopic 
width, is accomplished assuming that any spatial vari- 
ation of the irradiance can be neglected. This is valid 
because the laser spot diameter is assumed to be much 



larger than the region under study. The irradiance is 
thus, in effect, spatially constant over the infinite x-y 
plane by virtue of the PBC. 

The relatively large number of photons contained in a 
pulse, typically a few tens of thousands, ensures a uni- 
form spatial distribution of the energy at the surface. 
The temporal Gaussian distribution is simulated by the 
successive arrival of planes of photons spanning the en- 
tire surface of the supercell, the number of which being 
determined by the instantaneous irradiance; they are sep- 
arated in time by intervals ranging from At to typically 
10 x At. 



D. Absorption of light 

The absorption of light in silicon proceeds by the ex- 
citation of valence electrons via intraband and interband 
transitions, the frequency of which depends on several 
parameters, such as the density of free carriers and the 
wavelength and intensity of the laser pulse. The absorp- 
tion of photons of energy hv larger than the bandgap 
Eq can induce three types of optical transitions: (i) An 
interband transition following the absorption of a single 
photon by a valence-band electron; a bond is broken and 
the electron is promoted to the conduction band, leav- 
ing a hole behind. The newly created pair shares an 
energy ~ (hv — Eg), each carrier receiving an initial ki- 
netic energy determined by the set of selection rules, (ii) 
A nonlinear, two- (or more) photon, interband transi- 
tion, occurring when n > 2 photons are simultaneously 
absorbed. If nhv > <j>, where <j> = 4.85 eV is the silicon 
work function.cfj the electron can be photoemitted; the 
condition nhv < ef>, i.e., A > 255 nm when n = 1, thus 
ensures that no photoemission takes place, (iii) An in- 
traband transition, or free-carrier absorption, according 
to which the photon is absorbed by an electron already 
in the conduction band. 

Following their creation, the carriers, because of mo- 
mentum and energy conservation rules, occupy thin en- 
ergy shells and a small volume of k space. In a charac- 
teristic time T e h of typically a few tens of fs, the hole and 
electron subsystems achieve, through electron- hole scat- 
tering, a quasi-equilibrium state described by a Fermi- 
Dirac (FD) distribution at a common electronic temper- 
ature T e which is higher than the lattice temperature T. 
Thus, after ~ 100 fs, the carriers and the lattice con- 
stitute two decoupled subsystems, most of the energy 
remaining stored within the former (T e > T). 

If n represents the conduction band carrier density, the 
balance equation for e-h pairs can be written ascJ 

- ' - ' (2) 



dt 



V ■ J = G + R 



where J is the current density, G the carrier generation 
rate resulting from the absorption of photons, and R the 
net recombination rate given by 

,3 



R = -7?r 5 + 5{T e )n 



(3) 
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with 7 the Auger recombination coefficient and S(T e ) the 
impact ionization coefficient. The carrier generation rate 
G is the sum of one- and two-photon processes, G = 
G\ + G2, and can be written as 



G = 



(1-T)al(z,t) (1-T) 2 l31 2 (z,t) 



h v 



2hv 



(4) 



(1— T)I(z, t) being the intensity of the laser at depth z be- 
low the surface, T the reflection coefficient at the surface, 
and a and (3 the one-photon and two-photon interband 
absorption coefficients, respectively. The absorption pro- 
ceeds according to the Beer-Lambert law 



dl 

— = -(a + <dn)I-l3I 2 
az 



(5) 



where represents the free-carrier absorption cross- 
section. 

It is instructive to compare the two-photon to one- 
photon carrier generation rates, that is the ratio 



^2 = (i-r)/(2,i)A 

Gi 2a 



(0) 



If we choose A = 308 nm, each photon carries an en- 
ergy hv — 4.03 eV, laarer than the bandgap of silicon at 
300 K, E G ~ 1.12 eV.B Fog-Si at 300 K, V = T ~ (L59JB 
a = a ~ 1.5 X 10 6 cm" 1 ^ and ~ 40 cm/GW.N For 
an irradiance as high as 10 3 GW/cm 2 , a value at least 
an order of magnitude higher than the ones considered 
in this work, Eg- (||) gives G%jG\ < 0.01. Moreover, it 
has been shownE3 that free-carrier absorption is negligible 
at short wavelengths, typically below 1000 nm. There- 
fore, the one-photon interband transition is the domi- 
nant mechanism for absorption in the present case and 
all other processes are ignored; as noted above, this al- 
lows the photoemission of electrons to be neglected. 

In practice, upon arrival at the surface, the x-y position 
of a photon is determined at random following a uniform 
probability distribution and reflected, again at random, 
according toc3 



r(T) = r„ir fi xr, 



(7) 



whe re Tr = 4 x 10~ 5 . Upon melting of the surface (see 
Sec. [II A ), thcj-eflectivity is assumed to be that of /-Si, 
i.e., Ti — 0.73J13 If not reflected, the photon is absorbed 
in agreement with Eq. (^J), neglecting two-photon and 
free-carrier absorption, at depthEll 



z = —a 1 ln(£) , 

where £ is a number generated from a unifora 
tion; a is the absorption coefficient, given bjE 



where 0i 



ct(T) = ag exp 
300 K and 6 2 = 



T — 0i 

©2 

4680 K. 



(8) 

distribu- 
(9) 



The silicon atom closest to the absorption site, and 
which possesses at least one valence electron, is then ex- 
cited, i.e., its number of valence electrons is decremented 
by one. An e-h pair is created with an initial position 
given by the (x, y, z) coordinates of the excited atom. 
The relatively large absorption coefficient at A = 308 nm, 
combined with a depth of 60 nm for the simulation box, 
ensures that more than 99% of the energy is absorbed 
within the supercell. 

The determination of the exact initial kinetic energy of 
each carrier requires the knowledge of the allowed opti- 
cal transitions. However, because very fast carrier-carrier 
scattering leads to a FD distribution at temperature T e 
within a few tens of fs, it is reasonable to assume the gas 
to be in a state of an instantaneous quasi-equilibrium. 
Further, because T e is relatively large, typically a few 
thousand degrees, the FD distribution is well approxi- 
mated by a MB distribution. Thus, as an electron is 
promoted to the conduction band, its initial kinetic en- 
ergy is calculated according to a MB distribution at the 
current temperature T e . The problem of determining the 
set of optical transitions as well as the carrier density 
and temperature-dependent chemical potential, required 
to characterize a FD distribution, can thus be avoided. 
Using a Box-Muller transformation,^ a MB distribution 
at T e is obtained by calculating the electron initial kinetic 
energy, E^ i , according to 



Ei- 



-knT P 



;in(C i )cos 2 (2 7 re i ) , 



(10) 



»=i 



where Q and £j are random numbers taken from a uni- 
form distribution. The hole thus receives an initial ki- 
netic energy equal to (hv — Eq — E^)- The relationship 
used for Eg (in eV) as a function of the lattice temper- 
ature T and the carrier concentration n isE3E3 



E G 



kT 



(T + A) 



(11) 



where E G = 1.17 eV, k = 4.730 x 10~ 4 eV/K, A = 636 K 
and <; = 1.5 x 10 -8 eV cm. 



E. Relaxation processes 

A number of microscopic mechanisms allow the sys- 
tem to restore equilibrium following irradiation by a laser 
pulse; these are carrier-phonon scattering, carrier diffu- 
sion, impact ionization and Auger recombination. 

Impact ionization is the process through which an e-h 
pair is created when an energetic electron collides with 
a valence electron: the former transfers part of its ki- 
netic energy to the latter which is promoted to the con- 
duction band. Because impact ionization most likely 
originates from the generation of highly- energetic carri- 
ers through two-phpion interband transitions (negligible 
here; see Sec. IID)J13 it can be ignored. 
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In silicon, the dominant recombination mechanism is 
Augen£3 for which the characteristic recombination time, 
due to. screening effects at high carrier densities, is > 
6 psO Because the carrier density gradient, which de- 
velops at the surface following laser irradiation, is taken 
into account, in our model, by assuming a unidirectional 
flow of el ectron s and holes along the z axis into the bulk 
(see Sec. HE 2 ), the latter have usually diffused out of 
the simulation box in a time less than 6 ps. As a result, 
Auger recombination can also be neglected. 

The important relaxation mechanisms, in the present 
context, are carrier-phonon scattering and carrier diffu- 
sion, which we explicitly take into account. 



If 

valui 



Carrier-phonon scattering 



.carrier density n remains below the critical 
J n r ~ 10 22 cm" 3 , relaxation is thermalBj i.e., 



the carriers relax by transferring their kinetic energy to 
the lattice by the spontaneous emission of (mainly) op- 
tical phonons in a characteristic time t^q ^ 1 ps (also 
known as the carrier-phonon energy relaxation time te). 
The population of optical phonons thus increases for ap- 
proximately 1 ps, after which the latter relax by emitting 
acoustic phonons in a characteristic time tta ~ 10 ps. 
Eventually, the energy is redistributed among all the 
vibrational modes of the lattice and a state of quasi- 
cquilibrium is reached, characterized by a Bose-Einstein 
distribution at temperature T. The further evolution of 
the system can then be described through thermal pro- 
cesses only as T — > T e , but T is still higher than its 
initial value. 

There are several types of carrier-phonon scattering 
processes. It is possible to identify the most important 
ones, as well as their associated mean scattering rates. 
The dominant mechanisms for energetic electrons are the 
/ and g intervalley scattering transitions, which involve 
large changes in momentum and thus acoustic or optical 
phonons with wave vectors near the zone boundary.E3 

The probability P^t of observing a given intervalley 
transition, according to which a carrier in the valley i, 
with crystal momentum p^, scatters to a state with crys- 
tal momentum P/ in the valley /, is proportional to 
the square of the intervalley deformation potentiaJ T ..Dj f , 
which characterizes the strength of the scattering.E3 The 
probability can then be estimated as the weight associ- 
ated with Dfj . The probabilities thus obtained are listed 
in Table [| ancL-were calculated using the corresponding 
values for Z>;/.E3 As one can see, nearly 70% of phonons 
that scatter with electrons are optical. For holes, how- 
ever, the main_mechanism is optical deformation poten- 
tial scatteringjEj and we thus assume tha-holes to scatter 
only with optical phonons of 62.6 meV.E3 Moreover, the 
energy of the carrier determines whether scattering in- 
volves the absorption or the emission of a phonon. In 
Si, phonon emission is largely dominant for carrier en- 
ergies exceeding E s ~ 50 meV, while phonon absorption 



dominates below E S BB 

To determine the rate at which the energy is trans- 
ferred from the electronic to the ionic degrees of freedom, 
the mean carrier-phonon scattering rate as a function of 
carrier energy, 1/t~,, is also needed. This quantity was 
calculated in Si at room temperature usins-the full band- 
structure of the material by Fischetti et a/.E3 A fit to their 
data was embedded in our computer model; variations 
with temperature were ignored. 

Because of screening effects, however, the observed 
carrier-phonon scattering rate decreases for carrier den- 
sities n > 10 21 cm -3 , as recently observed by Sjodin et 
aZ. pJpfcheir results confirmed theoretical predictions by 
Yoffa.E3 Screening effects were thus implemented in our 
model as suggested by Yoffa, with the new carrier-phonon 
scattering rate given by 



\n) 



l/r° 

/ CP 



(i + (^;) 2 ) 



(12) 



where n sc ~ 1.2 x 10 21 cm~ 3 is the critical carrier density 
to observe significant screening effects; r° is the value 
of T cp without screening. Finally, the probability of an 
electron suffering a collision with a phonon was assumed 
to be Drude-like: the probability, P s , that a carrier has 
scattered with a phonon after a time t, where t is the 
elapsed time since the last scattering event, is thus 



P s (t) = 1 - exp(-t/r cp ) 



(13) 



In practice, we proceed as follows: after being gen- 
erated, a carrier is given a n ini tial position and kinetic 
energy as described in Sec. II D and its motion proceeds 
until a scattering event takes place. At each time step, 
the time elapsed since the last scattering event is de- 
termined, r cp is computed and the probability P s is de- 
termined from Eq. (13|). If a scattering event occurs, a 
phonon is either emitted or absorbed depending if the 
carrier kinetic energy is greater or lower than E s , respec- 
tively. 

For an electron scattering with a phonon, the inter- 
valley process is determined at random according to the 
probabilities of Table |. A quantum of energy equal to 
Tiwq is then given to (phonon emission) or removed from 
(absorption) the lattice according to a spatial Gaussian 
distribution in a radius of 5 A from the carrier, the lat- 
ter value being an estimate of the carrier wave packet 
width based on the uncertainty principle. Though a £U 
nite duration is associated with each scattering event 
the latter are assumed to be instantaneous for simplicity. 
Finally, if a phonon scattering event has occurred, the 
carrier energy is updated. 



2. Carrier diffusion 

Carriers diffuse into the bulk as a result of the den- 
sity gradient, Vn, which is itself a consequence of the 
exponentially-decreasing absorption with depth. The 
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variations with n of the ambipolar diffusion coefficient 
D and carrier mobilities, which characterize the dynam- 
ics of carriers in semiconductors, have been the subject 
of controversy; the available experimental data at den- 
sities exceeding 10 19 ifm,~ 3 are very scarce and concerns 
mostly germanium.E3HEj In fact, for n > 10 19 cm~ 3 , it 
is not clear whether the ambipolar diffusion coefficient 
decreases due to additional carrier-carrier scattering, as 
proposed by Fle^jherjS or increases due to many-body 
quantum effec 
been proposec 



the same time, various models have 
to explain the influence of carrier- 
carrier scattering, as well as ajjantum effects, on D. 

Recent experimental dataE3'E§ for moderate (10 15 — 
10 17 cm -3 ) carrier densities tend to support Fletcher's 
model. In view of this, and the lack of experimental evi- 
dence on many-body quantum effects for n > 10 19 cm~ 3 , 
we adopt Fletcher's model for D implemcntcdriji our code 
using the expression suggested by Berz et alE3 

Thus, each e-h pair is given an ambipolar diffusion 
coefficient which depends upon the local carrier density 
and local temperature, and which is updated at every 
time step, during which a carrier travels a distance Az — 
%/ DAt away from the surface. 



III. RESULTS AND DISCUSSION 

All simulations were carried out assuming a laser pulse 
duration tl = 10 ps and photon energy hv = 4.03 eV 
(A = 308 nm). The pulse intensity has a temporal Gaus- 
sian profile, but is spatially uniform. The fluence, F, 
was varied between 0.01 and 0.75 J/cm 2 ; corresponding 
intensitie s are in the range 1-75 GW/cm 2 . 

In Sec. Ill A, we present our results for fluences of 0.01 
to 0.20 J/cmT No matter removal is observed; a breath- 
ing mode of the surface reveals the propagation of pres- 
sure waves in the solid. At and above F t h = 0.25 J/cm 2 , 
large pieces of molten material are expelled from the sur- 
face; the mechanisms by which they are ejected are ex- 
amined in Sec. IIIB. 



A. Below the threshold energy for ablation 

Fig. ^| shows the carrier temperature T e during the first 
two picoseconds for a fluence F = 0.05 J/cm 2 ; T e is ob- 
tained by summing the kinetic energies of the electrons 
and holes in the conduction and valence bands, respec- 
tively. The time evolution of the carrier density n and 
the surface, T s , and bulk, Tb, values of the lattice tem- 
perature are shown in Fig. [| for the same fluence; here, 
T s is obtained by averaging over the first 8 = 1/a *~ 7 nm 
below the surface, where S is the optical skin depth, while 
Tb represents an average over the rest of the supercell. 
The pulse starts at t — and the initial temperature of 
the target is ~ 300 K. Upon arrival of the first photons, 
T e rises rapidly to ~ 11000 K but n remains low: the 
pulse provides sufficient energy to create highly-energetic 



e-h pairs while not significantly increasing the number of 
carriers. The newly-generated hot electrons and holes 
promptly begin to thermalize with the lattice: within 
~ 500 fs to 1 ps, T e drops to ~ 1500 K; this indicates 
a very fast initial cooling rate of the carriers through 
phonon scattering. For delays > 1 ps, T e decreases at a 
slower rate, i.e., the rate at which energy is being trans- 
ferred from the electronic to the ionic degrees of freedom 
diminishes significantly. The carrier-phonon energy re- 
laxation time, te, is thus estimated to be a few hundred 
fs. These observations are consistent with recent prob- 
ing o-Lultrafast carrier dynamics in silicon by Goldman 
et alE3 who reported a value of te ~ 1 ps. Equilibrium 
with the lattice is finally achieved after ~ 10 ps with 
T e — ► 300 K. 

The carrier density n in Fig. |3|a reaches a peak at t ~ 
7 ps, thus well after T e has passed its maximum; this is 
because the e-h plasma loses energy to the lattice faster 
than it gains energy from the pulse. The decrease of n 
for t > 7 ps is due to diffusion of the e-h pairs away from 
the surface; moreover, the carrier density remains well 
below the critical value, n c ~ 10 22 cm~ 3 , ensuring that 
the TAM is appropriate. 

The pulse generates a high number of electrons and 
holes which, in turn, heat up the lattice; this is seen in 
Fig. ||b. For t > tl, Tb reaches a plateau at ~ 1250 K; 
this value is below the melting temperature of c-Si, 
T m = 1685 K.c3 The surface temperature T s , on the 
other hand, continues to increase beyond tl, approaching 
Tb asymptotically. Most importantly, T s < Tb for the 
whole duration of the simulation. This clearly indicates 
that the maximum lattice temperature is located beneath 
the surface, a direct consequence of carriers diffusing into 
the bulk as shown below. 

The fact that the maximum temperature is found be- 
low the surface can be better appreciated from Fig. ^, 
where the lattice temperature T is plotted as a function 
of depth z at t — 50 ps for a fluence F = 0.20 J/cm 2 ; 
z = corresponds to the surface. The temperature in- 
creases from a minimum of ~ 2000 K at the surface to a 
maximum of ~ 4000 K at z ~ 27 nm; the latter is above 
the boiling point T{, = 2753 K, butr-hclow the critical 
temperature, T c — 5193 K, of silicon.ca For z > 27 nm, 
T decreases. Thus, the maximum temperature is not oc- 
curring at the surface but, instead, at z ~ 30 nm. These 
observations can be explained as follows: the pulse in- 
tensity decreases exponentially with depth and, conse- 
quently, more e-h pairs are created near the surface (see 
Sec. p D| ). However, the time interval between two con- 
secutive carrier-phonon scattering events being typically 
~ 10 fs to a few ps, a carrier will diffuse, on average, a few 
nm away from the surface between two scattering events. 
Thus, a large proportion of the energy is not released at 
the surface but, rather, deeper in the bulk. 

Clearly, the temperature in Fig. ||| does not decay ex- 
ponentionally with depth, as one might have expected; 
rather, a steep temperature gradient has developed at 
z < 30 nm. This is a clear indication of SSSH effects 
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which are usually assumed to, originate from a vaporiza- 
tion process at the surface.EJ However, it is important 
to note that, in our simulations, no vaporization is ob- 
served; this contrasts, for example, with the results of 
Zhigilei et al. who identified evaporation as the mecha- 
nism of paatter removal below the threshold energy for 
ablationEcj The absence of evaporation in the present 
simulations is not surprising if one computes the flux (in 
atoms m~ 2 s-. 1 ) of atoms that desorb from the surface, 
$ e , given byO 



$e = 



y/2TTMk B T 



(14) 



where M is the mass of a desorbed atom and P e is the 
pressure of the vapor in equilibrium with the solid at a 
temperature T e s . For F > 0.15-J/cm 2 , T s ~ 2000 K 
(see Fig. g), P e ~ 2.11 Pa,Llf and we find $ e ~ 
2.3 x 10 22 atoms m~ 2 s -1 . The surface of the super- 
cell has an area 3x3 nm 2 = 9 x 10~ 18 m 2 and it is thus 
estimated that an atom desorbs every ~ 5 fis, which is 
much longer than the duration of our simulations. Thus, 
the temperature profile in Fig. |] cannot be explained by 
an evaporation process; rather, it is the result of carriers 
diffusing into the bulk. 

Fig. |5j shows the surface temperature T s as a function 
of fluence. The threshold for melting, F m , is thus esti- 
mated to be in the range 0.05-0.1 GW/cm 2 since we must 
have T s > T m to observe melting. A more precise assess- 
ment of F m is obtained by computing the coordination 
number averaged over the first few atomic layers below 
the surface, po, as a function of time for various fluences. 
The results for F = 0.09 J/cm 2 are depicted in Fig. 0. 
A sudden change indicating the onset of melting is ob- 
served at t ~ 75 ps: from ~ 3.88 (a value slightly below 
that of bulk c-Si, a consequence of including in the cal- 
culation undercoordinated surface atoms) , rpa approaches 
the coordination number of Z-Si, i.e., ~ 6.113 Because no 
significant changes in the surface coordination number 
are observed at lower fluences, the predicted threshold is 
thus F m ~ 0.09 J/cm 2 , melting taking place on a char- 
acteristic time of a hundred ps, in good agreement with 
experimental data.oEjO 

The mechanisms by which the system restores inter- 
nal equilibrium below F t h are revealed in Fig. where 
the surface position, i.e. the z coordinate of the top- 
most atomic layer, z top , is plotted with respect to time 
for F = 0.10 J/cm 2 ; z top = corresponds to the ini- 
tial surface position. Also plotted is the bulk pressure 
Pb, averaged over the entire supercell. As one can see, 
oscillations with an amplitude of ~ 1 nm and a period 
~ 30 ps are observed at the surface. This breathing mode 
is due to the propagation of pressure waves in the solid; 
the latter represent the elastic response of the material 
to the intense heating by the pulse. Upon inspection of 
Fig. (?| Pb and z top are clearly out of phase: because 
silicon expands when heated, Zt op initially increases. At 
t ~ 20 ps, z top is maximum and P B is minimum (Pg < 0, 
tension); this is indicative of tensile stresses. In order to 



restore internal equilibrium, the material subsequently 
enters a compressive phase and contracts, i.e., z top de- 
creases. At t ~ 35 ps, Zt op is minimum and Pb is max- 
imum (Pb > 0, compression); the solid is then under 
compressive stresses. It is also apparent from Fig. ^ that 
melting occurs at the surface: the average value of Zt op 
decreases with time, an indication that the material is 
more compact— Indeed, molten silicon is denser than crys- 
talline siliconJij 

The surface velocity, Vt op , and acceleration, cit op , com- 
puted from the data in Fig. 0, are found to be ~ 100 m/s 
and ~ 10 13 m/s 2 , respectively. These values are con- 
firmed by noting that, upo*^ heating, a surface undergoes 
an expansion Az given by 



Az=(l- T)Fa T /pC 



(15) 



where T is the reflectivity, ay is the linear thermal ex- 
pansion coefficient, p is the density and C is the sr. 
cific heat. Assuming T ~ 0.59, a T ~ 2-S.x 10~ 6 K" 
p ~ 2.33 g/cm 3 , C ~ 0.713 J/gK- 1 ^ and a fluence 
F ~ 0.10 J/cm" 2 at 300 K, Eq. @ yields an oscillation 
amplitude Az ~ 0.6 nm, a value consistent with the re- 
sults of Fig. [?]. With a pulse duration tl — 10 ps, the 
velocity and acceleration of the surface become Vt op ~ 
Az/tl ~ 60 m/s and a top ~ Az/r| ~ 10 13 m/s 2 , 
thus in good agreement with the above results. In ad- 
dition, Fig. U gives the maximum value of w top and 
ato P as a function of laser fluence. A maximum is ob- 
served at F = 0.15 J/cm 2 . The subsequent decrease at 
F = 0.20 J/cm 2 might result from a stronger damping of 
the oscillations due to a thicker molten layer at the sur- 
face. Laser-induced oscillations at surfaces have found 
several applications. In dry laser cleaning, for instance, 
the removal of contaminants, i.e., adsorbed micron-sized 
particulates, is driven by momentum transfer: though of 
small amplitude (~ 1 nm), the surface oscillations take 
place on a very short time scale and thus exhibit very 
large accelerations. 

Finally, Fig. g shows the pressure P as a function of 
depth z at t = and t = 10 ps; the fluence is F — 
0.10 J/cm 2 . As expected, the pressure is zero throughout 
the material before laser irradiation. By the end of the 
laser pulse, a pressure gradi ent has built up at 35 < z < 
40 nm; we show in Sec. [II B that the latter is responsible 
for matter ejection at fluences F > 0.25 J/cm 2 . Below 
the threshold energy for ablation, however, the pressure 
gradient is not high enough to cause ejection of matter, 
and the system relaxes through the generation of pressure 
waves. 



B. Above the threshold energy for ablation 

The threshold energy for ablation, Fth, is defined in 
this work as the minimum fluence for the onset of mat- 



ter removal. No evaporation is observed (see Sec. Ill A ) 
rather, a collective ejection process is obtained at all flu 
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ences above F t h- Here, 
to a laser intensity I t h 



F t h = 0.25 J/cm 2 , corresponding 
= 25 GW/cm 2 . 



Fig. [L0| shows the surface position, z top , as a function 
of time for a fluence F — 0.35 J/cm 2 . As is the case 
for F < Fth, the material initially expands upon heat- 
ing. However this does not result, this time, in a sur- 
face breathing mode. Instead, a dramatic change occurs 
at t ~ 30 ps: from a value of ~ 7 nm, z top drops to 
~ —35 nm within a few ps. This indicates that ablation 
has taken place over a depth h ~ 35 nm, a value con- 
sistent with experiment. Ej No significant change of h is 
observed at higher fluences. 

A snapshot of the ablation process at t = 33 ps for 
the same fluence is given in Fig. [□]. The photons were 
incident from the top. The atoms that have absorbed at 
least one photon are colored in dark gray, others are in 
light gray. Evidently, most photons have been absorbed 
near the surface. 

Two features should be noted: (i) Ablation is a bulk 
phenomenon, i.e., the ejection of matter originates in the 
bulk and not at the surface as one might have expected. 
More specifically, ablation results in the expulsion of a 
large piece — a few tens of nm thick — of molten ma- 
terial from the surface with an initial axial velocity of 
~ lOOOm/s. (ii) As opposed to a disordered network in 
the bulk, a few nm thick crystalline layer is visible at the 
surface. The thickness of the latter corresponds to the 
average distance, about 5 nm, carriers created at the sur- 
face travel before suffering a first collision. For z > 5 nm, 
the temperature rises rapidly and, as shown below, the 
thickness of the layer removed during ablation depends 
upon the location of the pressure build-up in the bulk. 

Insight into the mechanisms responsible for ablation is 
found in Fig. [l2] showing the pressure, P, and the lat- 
tice temperature, T, as a function of depth for F = 
0.35 J/cm 2 . A steep temperature gradient has devel- 
oped at the surface by t = 10 ps, Fig. |T^a: from a value 
of ~ 2000 K at z — 0, T rises rapidly to reach a plateau 
slightly above the c ritica l temperature, T c , for z ~ 30 nm. 
As shown in Sec. Ill A, the occurrence of a maximum 



temperature below the surface is not a consequence of 
evaporation from the surface, nor does it result from a 
thermal diffusion process carried by phonons. Rather, 
the energy is transported into the bulk by the electronic 
degrees of freedom; this mechanism is, of course, absent 
in organic solids B 

Intense heating by the pulse leads to the thermal con- 
finement of the laser-deposited energy and, as shown in 
Fig. fl2|a, the material is heated far beyond its boiling 
point by the end of the pulse. At the same time, the 
lattice cannot undergo sufficient thermal expansion and 
a strong pressure (> 0) in excess of 10 GPa develops at 
35 < z < 50 nm, Fig. [jjb. By t ~ 24 ps the latter has re- 
laxed, leaving important tensile stresses in the material; 
these are responsible for the rapid expansion observed for 
t > 10 ps in Fig. [l0|. The pressure profile at t — 24 ps in 
Fig. |l2| b also reveals a local minimum (in absolute value) 
at z ~ 35 nm. As a consequence, strong tensile forces of 



opposite directions drive the material apart locally and 
nucleation of voids occur, as depicted in Fig. [l3| The 
latter begins to form at t = 27 ps and ablation follows a 
few ps later as the mechanical strength of the material is 
exceeded. 

As seen above, the high compressive stress which de- 
velops within the absorbing volume is responsible for the 
subsequent tensile forces; thus, it is the dri ving force re- 
sponsible for matter removal above F t h- Fig. |14| shows the 
variation of the maximum compressive stress, P ma x, with 
laser fluence. Below the threshold for ablation, P max in- 
creases with F. Above Fth, however, JP m ax saturates. 

In a recent study of organic solidsjj Zhigilei et al. re- 
ported the observation of void nucleation in the stress 
confinement regime and attributed it to the build-up 
of a strong pressure gradient within the absorbing vol- 
ume. The condition of stress confinement is expressed 
by tl < ts ~ S/vs, where vg is the velocity of sound 
in the material. For c-Si, vs ~ 9 x 10 3 m/s and,c3 for 
a laser pulse at A = 308 nm, t$ ~ 1 ps. This value 
is < T£ by a factor of 10; conditions of stress confine- 
ment arc thus unlikely to develop in silicon following ir- 
radiation with 10 ps or longer pulses. The regime of 
thermal confinement, on the other hand, is expressed hv, 
tl < i~th ~ S 2 /Dt, where Dt is the thermal, diffusivity.Q 
For silicon at 300 K, D T ~ 0.86 cm 2 /s;ED this corre- 
sponds to r t h ~ 0.5 ps, a value which is also about an 
order of magnitude lower than tl- However, Dt drops 
to - 0.10 cm 2 /s as T reaches ~ 1500 K,0 and it follows 
that tl ~ T t h- this suggests that conditions of thermal 
confinement are likely to arise when silicon is heated be- 
yond T ~ 1500 K in a typical time of 10 ps or less. This 
phenomenon is responsible for the rapid heating of the 
material up to T ~ T c in Fig. |l2|a. However, one should 
note that the presence of a steep temperature gradient 
at the surface rules out the possibility of phase es 
sion mechanisms which require that dT/dz ~ 0J15 
Consequently, ablation is here the result of SSSH effects 
brought about by the thermal confinement of the energy. 

Thus, the present simulations suggest that picosecond 
laser ablation of silicon is driven by SSSH effects. How- 
ever, a more precise assessment of the matter removal 
mechanisms should involve the implementation of sur- 
face effects; these have been observed|iii a study of hot 
carrier dynamics at a Si(100) surfaceE3 Accounting for 
surface effects is currently under development. 



IV. CONCLUDING REMARKS 

A molecular-dynamics thermal annealing model 
(MADTAM) has been applied to the study of laser ab- 
lation of silicon with picosecond pulses on a 100 ps time 
scale. A detailed description of the thermal annealing 
model on a microscopic scale has been embedded into 
a molecular-dynamics scheme; this was accomplished by 
explicitly accounting for carrier-phonon scattering and 
carrier diffusion. 
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The model predicts a melting threshold F m = 
0.09 J/cm 2 and an ablation threshold F th — 0.25 J/cm 2 . 
Below Fth, oscillations of the surface reveal the propa- 
gation of pressure waves in the solid. Above F t h, im- 
portant subsurface superheating (SSSH) effects are re- 
sponsible for matter removal: as a result of the thermal 
confinement of the laser-deposited energy, the material is 
overheated to a temperature corresponding to the critical 
point of silicon and a pressure gradient builds up a few 
tens of nm beneath the surface. Ejection of matter is ini- 
tiated by the relaxation of the strong compressive stresses 
after a few tens of ps, large pieces of molten material be- 
ing expelled from the surface with high axial velocities. 

We are currently investigating the role of surface ef- 
fects. These are believed to be responsible for the for- 
mation of a surface space charge layer which is expected 
to influence the carrier dynamics and to compete with 



carrier diffusion. 
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TABLE I: Probabilities Pif associated with the intervalley 
processes. 

hup (meV) D if (10 8 eV/cm) P if 
Ts 62.6 (LO,TO)° L75 75 0.344 

gi 44.3 (TA) 6 1.18 6 0.156 

g 2 22.1 (LA) 6 1.18 b 0.156 

33 62.6 (LO) a 1.75 6 0.344 

Total 1.000 



a Ref. 
6 Ref. 
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FIG. 1: Generation of an e-h pair and the main subsequent 
relaxation mechanisms in MADTAM. (a) creation of an e- 
h pair upon absorption of a photon; (b) scattering of the 
electron or hole by a phonon; (c) unidirectional diffusion of 
the e-h pair into the bulk along the z axis. 
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FIG. 2: Carrier temperature T e as a function of time for a 
fiuence F = 0.05 J/cm 2 . 




FIG. 3: (a) Carrier density n and (b) surface T s and bulk 
Tb temperatures as a function of time for a fiuence F = 
0.05 J/cm 2 . 
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FIG. 4: Lattice temperature T as a function of depth at t 
50 ps. The fluence is F — 0.20 J/cm 2 . 
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FIG. 5: Surface temperature T s as a function of laser fluence. 




FIG. 6: Average surface coordination number po as a function 
of time for a fluence F = 0.09 J/cm 2 . 
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FIG. 7: Surface position zt op and bulk pressure Pb as a func- 
tion of time for a fiuence F = 0.10 J/cm 2 . 
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FIG. 8: Maximum velocity and maximum acceleration of the 
surface for different fiuences. 
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FIG. 9: Pressure P as a function of depth at different times. 
The fluence is F = 0.10 J/cm 2 . 
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FIG. 10: Surface position z top as a function of time. The 
fluence is F — 0.35 J/cm 2 . 



FIG. 11: Snapshot of the target at t — 33 ps illustrating mat- 
ter removal. Atoms that have absorbed a photon are depicted 
in dark gray; others are in light gray. The lateral and verti- 
cal dimensions are ~ 3 nm and ~ 40 nm, respectively. The 
fluence is F = 0.35 J/cm 2 . 
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FIG. 12: (a) Lattice temperature T as a function of depth at 
t = 10 ps and (b) pressure P at different times as a function 
of depth. The fluence is F = 0.35 J/cm 2 . 
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FIG. 13: Snapshots of the subsurface region at a depth 30 < 
z < 40 nm. The surface (not shown here) is above. Atoms 
that have absorbed a photon are depicted in dark gray; others 
are in light gray. The strong tensile stresses are responsible 
for void nucleation beneath the surface. The fluence is F = 
0.35 J/cm 2 . 
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FIG. 14: Maximum compressive stress Pmax below the surface 
as a function of fluence. 



